Dietary wheat and reduced methane yield are linked to rumen microbiome changes in dairy cows

Fermentation of pasture grasses and grains in the rumen of dairy cows and other ruminants produces methane as a by-product, wasting energy and contributing to the atmospheric load of greenhouse gasses. Many feeding trials in farmed ruminants have tested the impact of dietary components on feed efficiency, productivity and methane yield (MeY). Such diets remodel the rumen microbiome, altering bacterial, archaeal, fungal and protozoan populations, with an altered fermentation outcome. In dairy cows, some dietary grains can reduce enteric methane production. This is especially true of wheat, in comparison to corn or barley. Using a feeding trial of cows fed rolled wheat, corn or barley grain, in combination with hay and canola, we identified wheat-associated changes in the ruminal microbiome. Ruminal methane production, pH and VFA concentration data together with 16S rRNA gene amplicon sequences were used to compare ruminal bacterial and archaeal populations across diets. Differential abundance analysis of clustered sequences (OTU) identified members of the bacterial families Lachnospiraceae, Acidaminococcaceae, Eubacteriaceae, Prevotellaceae, Selenomonadaceae, Anaerovoracaceae and Fibrobacteraceae having a strong preference for growth in wheat-fed cows. Within the methanogenic archaea, (at >99% 16S rRNA sequence identity) the growth of Methanobrevibacter millerae was favoured by the non-wheat diets, while Methanobrevibacter olleyae was unaffected. From the wheat-preferring bacteria, correlation analysis found OTU strongly linked to reduced MeY, reduced pH and raised propionic acid levels. OTU from the genera Shuttleworthia and Prevotella_7 and especially Selenomonadaceae had high anti-methane correlations. An OTU likely representing (100% sequence identity) the fumarate-reducing, hydrogen-utilising, rumen bacterium Mitsuokella jalaludinii, had an especially high negative correlation coefficient (-0.83) versus MeY and moderate correlation (-0.6) with rumen pH, strongly suggesting much of the MeY suppression is due to reduced hydrogen availablity. Other OTU, representing as yet unknown species from the Selenomonadaceae family and the genera Prevotella_7, Fibrobacter and Syntrophococcus also had high to moderate negative MeY correlations, but low correlation with pH. These latter likely represent bacterial species able to reduce MeY without causing greater ruminal acidity, making them excellent candidates, provided they can be isolated, for development as anti-methane probiotics.


Introduction
While farmed ruminants such as cattle can feed on various forage species and grains, they do so at a cost to the environment. A major concern for global warming is the greenhouse gas methane emitted by ruminants such as sheep, goats and cattle [1,2]. This gas is a natural product of the fermentation process occurring in the rumen of animals digesting feeds such as grasses and grains.
Many dietary interventions have been tested in farmed ruminants, such as dairy cows, aiming to reduce methane emissions or improve feed efficiency. Different pasture and forage types have been evaluated along with a variety of edible materials such as various grains [3], plant extracts and chemical additives such as 3-nitrooxypropanol [4][5][6]. The benefits and drawbacks of a variety of dietary additives including wheat and Asparagopsis seaweed have been recently reviewed [7].
An important theme seen in many of these studies has been the redirection of metabolic hydrogen away from production of methane and more toward propionate [8] or other reduction reactions. The type of grain fed may influence methane production from the rumen, but it remains unclear what impact diets have on this microbiome or how the effect may be caused. Is it due to starch degradation rate for example, or some unidentified factor in one grain or another? Perhaps also the outcome of feeding different diets may depend on the pre-existing state of the ruminal microbial population [9]. Different forages and grains have been compared and exotic additives tested in an effort to manipulate the fermentation process occurring in the rumen [6]. Several feeding experiments have included an investigation of the impact of dietary changes on the rumen microbiome [10]. Comparisons have been made using the commonly available grains such as corn, wheat, oats and barley, resulting in some promising, but sometimes contradictory, reductions in methane production. For example, in sheep, a diet containing wheat was found to generate more methane than a diet including corn [11], whereas in dairy cows the opposite is true. Moate et al. [12] found that inclusion of wheat in the diet reduced methane production significantly when compared to dietary inclusion of corn or barley. Here we used 16S rRNA sequencing to study the rumen microbiome of dairy cows fed, over a 3-week period, with either barley, corn or wheat. Understanding the rumen microbiome is key to understanding rumen fermentation and the changes it undergoes in response to diet. We hypothesise that dairy cows fed on diets supplemented with different grains would have different rumen microbiome populations and that such microbial differences would explain the lower methane production seen with the diet supplemented with wheat compared to barley or corn.

Ethics approval and consent to participate
All cows were maintained in the research herd at the Agriculture Victoria Research Ellinbank Centre, 1301 Hazeldean Road, Ellinbank, Victoria, Australia. All experiments were conducted in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes (National Health and Medical Research Council 2013). Approval to proceed was obtained from the DJPR Agricultural Research and Extension Animal Ethics Committee. No anaesthesia, euthanasia, or any kind of animal sacrifice was part of the study.

Diet
The rumen samples collected in this study are from the experiment described previously on ruminal methane [12]. Cows were fed on hay, canola meal and grains as described [12]. Essentially the cows were fed 55% hay + canola meal + minerals and supplemented with 45% rolled grain (either corn, barley or wheat). The timeline of the feeding trial is summarised in S1A Fig. Rumen samples for microbiome analysis were collected in parallel to those used for chemical analysis and stored at -80C until processing for DNA sequencing.
Thirty-two lactating, multiparous Holstein-Friesian cows, including 12 fistulated cows, were assigned to one of 4 blocks, each block consisting of 5 non-fistulated and 3 fistulated cows. Each block was balanced for cow age, weight and days in milk (DIM) according to Baird [13] then randomly assigned to one of four dietary treatments: Each day at 0700hr and 1600hr, each of the 32 cows was fed in its individual feed stall provided with its specified feed and a supply of water.
Cold-pressed canola meal was included in these diets to ensure that the diets met the cows dietary requirements for protein and that the diet also supplied an amount of fat (i.e. > 500g fat/cow/day) equivalent to at least 50% of the expected milk fat yield of the cows.
The first week of the experiment constituted a covariate period during which all cows were fed a common diet comprising 6 kg DM/cow/day of rolled wheat and 16 kg DM/cow/day of alfalfa hay. The second week of the experiment constituted a transitioning period during which cows were transitioned onto their experimental diets. The third week of the experiment was an adjustment period during which the cows were fed their full experimental diets. The fourth and fifth (last) week of the experiment was an intensive period of measurements. Samples of ruminal contents were collected for microbiome analysis on the last day of week 5. Methane measurements were taken over the last 6 days of the trial.
Note that cows cn2319, wt9534 and wt9543 were treated with antibiotics (CepravinLC, Yodimaspen and Mastalone respectively) for mastitis in the month leading up to the transition period, before beginning the experimental diet treatments.

Methane, pH and VFA measurements
Methane emissions were estimated using the SF 6 tracer technique as described [14]. Using this method, methane measurement requires the controlled release of a tracer gas, sulphur hexafluoride (SF6) into the rumen. Gelatin-coated capsules containing about 2.8g SF6 were inserted into the rumen on day 24 of the trial. An animal's methane emission was calculated from the measured gas-mixing ratio in eructated gases (SF6 tracer + CH4). Details can be found elsewhere [14]. Gas was collected using an evacuated canister connected to a suction tube near the cow's mouth and was sampled for 24 hours every day over the last 6 days of the trial. The average of the methane production and consumed feed (dry matter intake, DMI) measurements over the last 6 days of the experiment was used to calculate the average gm CH4 per kg DMI (MeY) for this analysis. See S1A Fig which shows a timeline of the diet treatment and methane measurements. The pH of ruminal samples was measured immediately after collection using a pH meter (Mettler-Toledo. Also measured were ammonia using the QuikChem method 12-107-06-2-A from Lachat Instruments, D-lactate and L-lactate determined by measurement of NADH produced by -L or -D lactate dehydrogenase [15] and acetic, propionic, iso-butyric, butyric, iso-valeric, valeric, hexanoic and heptanoic acids determined using capillary gas chromatography as described by Packer et al. [16].

Collection of ruminal content samples, genomic DNA extraction and 16S rRNA gene V4 amplicon sequencing
On day 35, the final day of the trial (21 days after the transition period), ruminal material was collected at approximately 11.30 AM from all cows using an oesophageal probe as described [4]. Rumen collection was by mouth using a suction tube. Two samples were collected, one frozen for bacterial microbiome analysis, one for chemical analysis. Samples (20ml) were kept frozen at -70˚C until analysed. Genomic DNA was extracted from a 200 μL aliquot from each sample using the PSP Spin Stool DNA Kit (STRATEC Molecular GmbH, D-13125 Berlin, Germany) according to the manufacturer's instructions.
Variable region 4 (V4) from the prokaryote 16S rRNA gene was amplified from the ruminal genomic DNA by polymerase chain reaction (PCR) using the F515/R806 primers (5'-GTGCCAGCMGCCGCGGTAA-3'/ 5' -GGACTACHVGGGTWTCTAAT-3') and associated published method [17], using Phusion DNA polymerase (ThermoFisher Scientific). The Illumina MiSeq system was used to generate DNA sequence data from the PCR fragments ( The no-DNA controls were samples with the same detergent, buffer, DNA polymerase enzyme and other required ingredients, but without ruminal DNA. Instead, the DNA was replaced with ultra-pure water. They were processed in parallel with the DNA from ruminal samples. The ruminal DNA samples and 2 no-DNA controls were all subjected to the same PCR reactions and DNA sequencing together as a batch.

Processing of 16S V4 amplicon DNA sequences
DNA sequence fastq files from MiSeq paired end (PE) sequencing were quality filtered using Trimmomatic software [18] using a sliding window method. Moving from one end of a sequence the 4-base window would cut the sequence if the average phred score of the 4 bases in the window was less than 15. If the remaining sequence was then less than 200 bases it was rejected. Quality score by base position plots of raw and filtered sequence data were generated by fastQC [19] and examined manually. Processing of fastq files, containing overlapping 5' or 3' end sequences, was carried out using basic unix utilities and the PANDAseq program [20] for assembly of paired-end sequences. Assembled sequence pairs were retained if they contained the 5' and 3' PCR primer sequences at their ends (using regular expressions based on the above sequences). Primer sequences were then removed from the sequence ends and sequences rejected if an ambiguous base (an N character) was present anywhere in the sequence. Unix utilities such as grep, gawk, sort and uniq were then used to aggregate and count replicate (100% identical) sequences, creating files of unique amplicon sequences of known abundance in fasta format.
All amplicon sequences found in the No-DNA controls were assumed to be contaminants and were removed from the ruminal sample sequences before any further analysis.

Rarefaction analysis of microbiome sequences
Groups of merged paired-end sequences were randomly selected using the unix utility shuf, in increasing multiples of 2000. The groups were clustered using Uparse [21] at 97% identity and the count of OTU (clusters) plotted (using R) versus the total sequences in each group.

Diversity and richness of microbiome sequences
Merged paired-end sequences were analysed for Shannon diversity index and Chao1 species richness using the R libraries vegan and fossil respectively. The values obtained were plotted for each ruminal sample.

OTU clustering
Filtered and merged 16S rRNA V4 paired-end sequences that had then been aggregated (at 100% identity) to create unique sequences, were further clustered to Operational Taxonomic Units (OTU) using the Swarm2 software [23] using parameters -f (fastidious) -d 1 (default). The abundance of each OTU was calculated using the abundance of each unique sequence and the total sequence count within each cluster (each cluster is known as a swarm). A swarm resembles a phylogenetic tree where each layer or branch of the tree differs from the one above by a single nucleotide and the central sequence or node is the centroid. The centroid is always the most abundant sequence of the cluster. For sequence comparisons (eg taxonomy) the centroid was used as the representative sequence and is referred to as the OTU. A fasta file of the 1260 most abundant OTU DNA sequences (swarms) is available at Figshare (Digital Object Identifier = 10.6084/m9.figshare.16689487) along with the MiSeq sequence files.

Statistical analysis and plotting
Matrices of samples by OTU or taxon counts were constructed using python version 2.7. All statistical analyses such as principal components (PCA), differential abundance, correlations (Spearman Rank Correlation) and the plotting of bar and volcano plots, PCA and correlograms were carried out using the R (ver 3.6) libraries Cairo, vegan, fossil, permute, lattice, latti-ceExtra, colorspace, RColorBrewer, calibrate, ggplot2, corrplot, ggrepel, ggfortify, grid, gridExtra, methods and/or edgeR. All R libraries except edgeR (for differential abundance analysis) were obtained from https://cran.r-project.org. The edgeR library was from https:// bioconductor.org/packages/release/bioc/html/edgeR.html. In the edgeR analysis, differential abundance was estimated using the quasi-likelihood F-test with the generalised linear model glmQLFit and glmQLFTest functions recommended by the edgeR authors [24,25] for differential expression calculations. The resulting fold changes (FC) are expressed as log 2 (FC). For PCA analysis (see below) the cows were assigned to MeY (g CH 4 /kg dry matter intake) groups as follows: low (MeY < 15), medium (15 < MeY < 25) or high (MeY > 25). These groups were based on the mean +/-a standard deviation. This is illustrated in S1B Fig.

Microbial families versus diet
The initial analysis of taxonomies present in the microbiome rRNA gene V4 amplicon sequences was based on families. Using the Silva 16S rRNA microbial taxonomy database, we found the sequence data from the combined sample collection (32 cows) represented 153 microbial families. At an abundance of at least 100 counts per million (cpm) in at least 5 cows (excluding the unassigned sequences) we found 46 microbial families: 2 archaeal families and 44 bacterial families. When considering the cows as diet-based groups of 8, no group was devoid of any of the 46 microbial families, despite relative abundances varying greatly. The 46 families and their abundances in cpm in each cow are shown in S2A Table. A core of 19 families (excluding unassigned), all bacterial, no archaeal, were identified in each of the 32 rumen samples at a minimum of 100 cpm abundance. These are indicated with "#" in S2A Table. The  Table). A basic visual interpretation shows firstly that, at the family level, the ruminal microbiomes of corn and barley-fed dairy cows were almost indistinguishable, being dominated, as expected, by the Prevotellaceae. In contrast to the corn and barley-fed cows, in cows fed the wheat based diet, a t-test shows there were substantial differences in abundance of some of the bacterial families (S2B Table). In cows fed the wheat diet, Acholeplasmataceae, Bacteroidales, Christensenellaceae, Defluviitaleaceae, F082, Hungateiclostridiaceae, Methanobacteriaceae, Muribaculaceae, Oscillospiraceae, p-251-o5, Prevotellaceae, Rikenellaceae, Ruminococcaceae, Spirochaetaceae, and UCG-011 were less abundant while Acidaminococcaceae, Lachnospiraceae, Veillonellaceae and Bifidobacteriaceae were proportionally greater in abundance in wheat-fed compared to corn or barley-fed cows.

OTU sequences
From the clustered 16S amplicon sequences we identified 7625 clusters or OTU. Of those OTU, 1260 had an abundance of at least 100cpm in at least 5 rumen samples. S3A Table contains a list of the 1260 OTU with their abundance in cpm in each cow's ruminal microbiome, plus the MeY for each cow. S3B Table contains a list of the same 1260 OTU with their likely taxonomy from the Silva 16s database.
The microbial populations were examined using principal component analysis (PCA), using both the distribution of microbial families and the distribution of OTU from the Swarm clustering. Fig 3A shows a PCA of the families and 3B the OTU. There was marginal clustering visible when basing the analysis on families, but a convincing grouping can be seen of OTU by diet and by MeY. The only two wheat-fed animals (wt9534 & wt9543) that did not cluster with wheat grouping in Fig 3B, and which had higher MeY, had both been treated with antibiotics in the lead up to the experiment. The latter outlier cows were ignored in subsequent analyses. Cow cn2319 was also treated with an antibiotic but did not behave as an outlier in the PCA. This signals that each diet results in a distinct microbial population able to influence MeY. In addition, the PCA confirms there was little if any difference between the two barley diets (single versus double rolled). In further analyses these were treated as a single group of 16 cows.
Of the 1260 OTU present at a minimum 100cpm in 5 cows, a core group of 36 OTU were present in the rumen of every cow (at a minimum 10cpm), although with abundance varying greatly. For example, OTU rs0037 ranged from 29cpm in cow cn2319 to 11854cpm in cow bs1226. Table 1 lists the 36 core OTU with their taxon path as determined by alignment of the DNA sequences versus the Silva 16S rRNA database, using Megablast. Again, Prevotellaceae and Lachnospiraceae family members made up the bulk of the core OTU and Methanobacteriaceae was the only archaeal representative.

Differential OTU abundance versus diet and methane yield
Differential abundance of OTU when contrasting diets, carried out using the edgeR library of R, identified many OTU detectable only when a specific diet was fed to the cows. Note that cows wt9543 and wt9534 were identified as outliers in the PCA and were ignored in this and subsequent analyses. The results from edgeR analysis can be visualised using a volcano plot  Table with p-value in 4b), creating a plot of MeY correlation versus dietinduced OTU abundance log 2 fold change (Fig 4), contrasting wheat with corn plus barley (ie non-wheat). The same figure but with OTU labels instead of Genera can be seen in S4B Fig. There were 395 OTU with mean abundance at least 128cpm and a significant positive or negative correlation with methane abundance (P < 0.01). We used these latter bacterial OTU, plus all the archaea-derived OTU, as input to the edgeR differential abundance analysis to identify wheat-linked OTU that are also linked to suppression of methane. Several OTU from the Acidaminococcus, Fibrobacter, Lachnoclostridium, Lachnospira, Mitsuokella, Prevotella_7, Pseudoramibacter, Selenomonas, Shuttleworthia and Syntrophococcus genera were significantly enhanced in the wheat-fed rumen microbiome and some correlate strongly with reduced methane, when compared to the non-wheat-fed rumen (Fig 4). Table 2 contains a list of the 25 OTU with the highest fold change in abundance in the wheat-fed rumen samples when comparing them to non-wheat diets. In the reverse direction, genera with high abundance in the corn or barley diets plus high positive MeY correlation include the Christensenellaceae, Flexilinea, Moryella, Prevotella, Prevotellaceae, Rikenellaceae, Sutterella plus unidentifiable members of the Rikenellaceae, Lachnospiraceae and Muribaculaceae families. Finally, the methanogens also have an important presence in these microbiomes. Several members of the Methanobacteriaceae and Methanomethylophilaceae families were increased in abundance in the nonwheat relative to the wheat-fed microbiome, while others of these families showed no difference between wheat and non-wheat diets.

Methanogens
The sequence clustering yielded 9 OTU from methanogenic archaea, belonging to the Methanobacteriaceae and Methanomethylophilaceae families. Two OTU dominated the archaeal population. Using the Genbank 16S rRNA database to obtain likely genus and species-level taxa, Methanobrevibacter millerae (rs0036, 99.6% sequence identity) was present more in the barley or corn-fed cows compared to wheat. In contrast, Methanobrevibacter olleyae (rs0037, 99.2% sequence identity) was present in all cows and its abundance had no significant correlation with MeY. The correlation analysis versus MeY gave rs0036 a moderate positive coefficient of 0.67(P = 0.001), while that of rs0037 is only 0.07 (P = 0.5). A one-tailed t-test confirms rs0036 abundance was greater in the rumen of corn or barley-fed cows than in cows fed wheat (see S5  Table).

OTU correlations versus MeY, pH and other fermentation products
Spearman Rank Correlation analysis of 1260 OTU (those with abundance at least 100cpm in at least 5 cows) in comparison to MeY levels identified 159 OTU with a positive MeY correlation greater than or equal to 0.7 and 31 OTU with a negative MeY correlation less than or equal to -0.7. The 83 OTU with an MeY correlation P-value less than 0.001, a logFC (differential abundance, wheat vs non-wheat) magnitude greater than 2.5 and an OTU abundance of 1000cpm in at least 5 cows were analysed for correlations with MeY, pH, NH4, D-lactate, L-lactate and acetic, propionic, iso-butyric, butyric, iso-valeric, valeric, hexanoic and heptanoic acids (S6 Table). The methanogen OTU rs0036 and rs0037 were included. The correlation coefficients can be seen plotted as a correlogram in Fig 5. Correlation results for the methanogens plus the 20 OTU with the most negative (ie most anti-methane) MeY correlation are listed in Table 3 (details in S4A and S4B Table), together with pH, lactate and those for acetic and propionic acids and their taxonomic assignments from the Silva database. Other factors analysed yielded only weak correlations.
The OTU rs0036, most likely from Methanobrevibacter millerae, had a strong correlation with MeY and pH plus a strong negative correlation with propionic acid. However, the OTU rs0037, most likely from Methanobrevibacter olleyae, had almost no correlation with MeY and no more than a weak correlation with any of the other factors measured. Of the bacteria, rs0090, most likely from Mitsuokella jalaludinii, had the strongest negative correlation, -0.83, with MeY. This OTU had a moderate negative correlation with pH together with a moderate positive correlation with lactate and a strong correlation with propionic acid. As the logFC suggests, OTU rs0090 is predominantly found in the wheat-fed rumen. Following on from rs0090, the OTU (rs0054, rs0111 and rs0215) from the Shuttleworthia genus and Selenomonadaceae family also had strong negative correlations with MeY, negative correlations with pH The 25 OTU found in the wheat-fed cows with the greatest change in abundance compared to corn or barley-fed cows, sorted by log 2 fold change. Sequence identity (ID) is based on alignment with 16S sequences in the Silva database.
https://doi.org/10.1371/journal.pone.0268157.t002 and positive correlations with lactate and propionic acid. Note that OTU rs0215, from the Selenomonadaceae family, of unknown genus or species, but having 98% sequence identity with Mitsuokella jalaludinii 16S rRNA V4, had a strong negative correlation with MeY but only a weak negative correlation with pH. OTU rs0109, also from Shuttleworthia has a strong correlation with MeY and is strongly associated with wheat, but is of lower abundance and so does not appear in this figure. More details can be seen in Table 3 and S4A and S4B Table. MeY itself had only a weak positive correlation with ruminal pH, but a moderate negative correlation with lactate and a strong negative correlation with propionic acid. Ruminal pH had strong negative correlations with acetic and propionic acids.

Antibiotics
S7 Table shows the results of z-score calculations to discover which OTU were depressed in cows administered antibiotics for treatment of mastitis prior to the trial ((Mastalone/wt9543, Yodimaspen/wt9534, CepravinLC/cn2319). Among those OTU with a z-score less than -1 are the OTU with the strongest MeY and pH negative correlations. In the Mastalone-treated, wheat-fed cow, wt9543, the 2 members of the Shuttleworthia genus (OTU rs0054 and rs0111) were very much depressed in abundance (1.8 and 2 standard deviations below the mean of the other wheat-fed cows). These 2 OTU have the strongest negative correlations with ruminal pH (-0.7 and -0.68 respectively). The ruminal pH itself in this cow had increased by more than 2 standard deviations from the mean of the pH in the other wheat-fed cows (from a mean of 5.89 to a pH of 7.14). The OTU rs0090, from Mitsuokella and having the strongest negative correlation with MeY (-0.83), was also depressed by Mastalone, by 1 standard deviation from the mean. In the Mastalone-treated cow MeY increased by over 10 standard deviations above the mean for the wheat-fed group.

Discussion
The first port of call when analysing the impact of a dietary treatment on the rumen microbiome is usually the archaeal methanogen population. While several different families of archaea were seen, mostly at low abundance, by far the most abundant were 2 members of the Methanobrevibacter genus. Based on a more than 99% 16S rRNA V4 sequence identity according to the Genbank 16S database, these were likely M. millerae and M. olleyae. These 2 species behaved differently. On average, the abundance of M. olleyae changed little from one diet to another, whereas M. millerae tended to be less abundant in the wheat-fed versus the corn or barley-fed cows. Coinciding with the difference in abundance, the OTU rs0036, from M. millerae, had a strong correlation with MeY and ruminal pH, consistent with previous studies [26] suggesting a link between low pH, M. millerae abundance and reduced methane yield. M. millerae abundance, MeY and pH all had very strong negative correlations with the level of propionic acid, but it is not clear whether propionic acid has a suppressing effect on MeY because of its acidity or because its production results in reduced hydrogen availability. This feeding trial provided us with fortuitous observations which may help elucidate possible associations between the use of antibiotics and effects on the rumen microbiome. The impact of antibiotic treatment on methane production and rumen microbiome makeup of the wheat-fed cows is important. While not carried out here as a deliberate experiment it suggests that key bacterial species, which would otherwise trigger a suppression of methanogen activity, are sensitive to such antibiotic treatment. The fact that the OTU with the greatest sensitivity to the antibiotics used also includes those with the highest correlation to MeY and ruminal pH means that in the future, with properly replicated experiments, antibiotics may constitute a valuable tool for dissecting the microbial population and identifying those key species.
The microbiome changes we observed when analysing bacterial and archaeal families in wheat-fed compared to corn or barley-fed cows is consistent with our hypothesis that microbial species changes in the wheat-fed rumen would explain the lowered MeY seen in those cows compared to the corn or barley-fed cows. Further analysis at the OTU level confirms the hypothesis. The principal component analysis (PCA), differential abundance and correlation analyses of 16S rRNA-based OTU abundances confirm that significant differences exist between the microbiomes of cows fed different grain-supplemented diets and that those microbial population changes are closely linked to MeY. This confirms our hypothesis regarding microbial species found predominantly in the rumen of wheat-fed cows being associated with a reduction in MeY. We found 86 OTU (representing different species or strains) with a very high fold change in abundance (more than 12 times higher in wheat than in corn or barley-fed cows). Of the latter, 25 were also strongly correlated to a decrease in MeY (Spearman Rank correlation below -0.7). High on the list of genera linked to the wheat-based diet and lower MeY, were genera such as the saccharolytic Shuttleworthia [27] and Prevotella_7 [28], associated with production of butyrate, acetate and lactate. Top of the list of wheat-preferring, anti-methane OTU, rs0090, likely representing the rumen bacterium Mitsuokella jalaludinii (100% 16S rRNA V4 sequence identity), from the Selenomonadaceae family, had the highest correlation (-0.83) with a reduction in MeY. M. jalaludinii is known to produce lactic and acetic acids [29], contributing to a lowering of pH, but more importantly it is able to reduce fumarate to succinate, which in turn can be reduced to propionate [30]. Interestingly, we found MeY had a stronger negative correlation to propionic acid than to either lactate, acetate or even to ruminal pH. The reduction reactions known to occur in M. jalaludinii would deprive the methanogenic archaea of hydrogen, a vital ingredient for methane production. Indeed, others have already shown the value of M. jalaludinii as a potential probiotic bacterial species for methane mitigation in an in vitro rumen fermentation system [31]. Other bacteria preferencing the wheat-supplemented diet, represented by OTU such as rs0215 from the Selenomonadaceae family and some also from the genera Prevotella_7, Fibrobacter and Syntrophococcus with high anti-methane correlation but low correlation to acidity, need more detailed molecular and functional analysis and likely represent important anti-methanogenic species with potential as probiotic treatments.
In conclusion we propose that wheat, as fed to dairy cows in our feeding trial, has as yet undefined properties which promote the growth of several different ruminal bacterial species seen only at very low abundance in cows fed other grains, such as corn or barley. The impact of these species is two-fold: firstly ruminal pH is lowered, likely restricting growth of one of the main methanogenic archaea (Methanobrevibacter millerae) and hence contributing to a lowering of MeY; secondly, and possibly of more importance in terms of MeY, growth of the rumen fumarate-reducing bacterium Mitsuokella jalaludinii is promoted, reducing the availability of hydrogen for the generation of methane and in the process, producing acetate, lactate and possibly propionate, contributing to ruminal acidity. Those OTU from species with high anti-methane correlation but low correlation to acidity, and hence a lowered acidosis risk, represent ideal candidate anti-methane probiotic species.   Table. Z-score analysis to determine which OTU were suppressed by antibiotics. Mastalone (cow wt9543), Yodimaspen (cow9534) or CepravinLC/cn2319). The list includes rs0111, rs0054, rs0090 and rs0215. These were in the top 10 OTU having the strongest anti-methane correlations. (DOCX)